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ABSTRACT 


The  purpose  of  this  paper  is  to  examine  the  computational  aspects  of  the 
minimum  fuel  continuous  low  thrust  orbit  transfer  problem  and  to  display- 
characteristic  numerical  features  introduced  by  various  physical  constraints. 
Minimum- fuel  orbit  transfer  by  low  thrust  is  typical  of  many  problems  in 
optimal  control  which  result  in  a  two  point  boundary  value  problem  which 
must  be  solved  by  some  iterative  numerical  procedure.  Two  techniques, 
Multiple  Substitution  Polynomials  (MSP)  and  Marquardt's  method,  are  shown 
to  be  applicable  to  this  task,  and  a  detailed  analysis  is  mace  of  the  behavior 
of  these  methods  in  the  context  of  the  low  thrust  problem.  A  variety  of  sub¬ 
problems  is  considered  with  parametric  variation  of  endpoints,  thrust-to- 
weight  ratio,  and  orbit  axial  orientation.  A  physical  barrier  is  found  which 
restricts  sample  points  in  certain  limiting  case  fixed  endpoint  transfers.  The 
existence  of  multiple  stationary  solutions  is  shown  for  the  case  of  intersecting 
orbits,  and  the  nearly  singular  behavior  in  that  region  is  investigated.  Numer¬ 
ical  results  for  several  transfers  are  found  to  compare  with  similar  results 
reported  elsewhere.  ( 
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I,  INTRODUCTION 


The  purpose  of  this  report  is  to  examine  computational  aspects  of  the 
minimum  fuel  continuous  low  thrust  orbit  transfer  problem  and  to  display- 
characteristic  numerical  features  introduced  by  various  physical  constraints. 
Computational  difficulties  in  finding  numerical  solutions  of  optimal  control 
problems  have  severely  limited  the  application  of  modern  control  theory.  The 
low  thrust  problem  is  typical  of  many  optimal  control  problems  for  which  the 
necessary  conditions  may  be  established,  and  yet  which  defy  analytic  solution 
(without  drastic  simplification)  or  straightforward  numerical  solution  because 
of  computational  difficulties. 

Orbit  transfer  by  low  thrust  has  received  considerable  attention  in  the 
past  decade  (Refs.  1  -  5),  and  numerical  solutions  have  been  obtained 
for  a  large  variety  of  subproblems  in  the  control  of  low-thrust  vehicles. 

However,  in  emphasizing  the  physical  characteristics  of  the  problem,  the 
computational  complexities  have  been  somewhat  overlooked.  There  has  been 
little  mention  of  the  behavior  of  various  numerical  techniques,  or  of  a  com¬ 
parison  of  techniques.  Very  often  the  schemes  employed  have  been  severely 
limited  in  general  application,  and,  in  most  cases,  the  extreme  sensitivity 
of  the  solution  to  the  choice  of  initial  conditions  has  dictated  knowing  the  form 
of  the  solution  a  priori.  The  two  point  boundary  value  problem  (TPBVP) 
arising  from  the  low  thrust  problem  has  been  solved  by  Newton- Raphson 
techniques  (Refs.  2,5),  quasilinearization  (Ref.  3),  and  dynamic  programming 
(Ref.  6),  although  the  computational  difficulties  encountered  have  usually  been 
formidable  for  all  but  the  simplest  of  orbit  ti*ansfer  conditions,  for  example, 
circle-to-circle  coplanar  transfers. 

In  this  report,  numerical  solutions  are  obtained  for  ellipse-to-ellipse 
transfers,  with  several  combinations  of  fixed  and  free  endpoints.  Of  particular 
interest  is  the  minimum  fuel  transfer  between  given  elliptical  orbits  in  which 
both  endpoints  are  free.  Optimal  solutions  are  obtained  numerically  for  transfers 
between  rotated  coplanar  elliptical  orbits,  and  the  existence  of  multiple 
extrema  is-  demonstrated.  Two  finite  dimensional  optimization  techniques, 
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Multiple  Substitution  Polynomials  (MSP)  and  Marquardt' s  method,  will  be 
examined  in  the  context  of  the  low  thrust  orbit  transfer  problem.  Attention 
will  be  focused  on  obtaining  a  successful  iterative  convergent  process  with 
emphasis  on  the  numerical  algorithm,  the  method  of  reduction  of  the  problem 
for  computational  purposes,  and  the  behavior  of  neighboring  solutions  for 
parametric  studies. 

The  problem  to  be  considered  is  (refer  to  Fig.  1)  to  determine  the 
optimal  thrust  orientrtion  history  to  transfer  a  constant  magnitude,  low-thrust 
vehicle  by  having  it  thrust  continuously  from  an  initial  orbit  to  some  desired 
terminal  orbit,  with  the  minimum  expenditure  of  fuel.  The  initial  and  terminal 
orbits  are  assumed  to  lie  in  a  common  plane  about  a  spherical,  drag-free 
earth.  Thus,  the  only  forces  acting  on  the  vehicle  are  its  thrust  and  its  own 
weight.  Plane  polar  coordinates  are  chosen  for  the  geometry  of  the  problem, 
with  the  thrust  orientation  angle  measured  clockwise  from  the  circumferential 
direction.  The  assumption  of  a  constant  mass  flow  rate  transforms  the  mini¬ 
mum  fuel  problem  to  a  minimum  time  problem,  with  fuel  expenditure  linearly 
proportional  to  total  transfer  time. 

The  orbits  chosen  have  perigee/apogee  altitude  values  of  55/110  and 
80/180  n  mi.  Initially,  the  orbits  are  taken  as  co-apsidal  (Aco=  0  deg),  but 
results  are  obtained  for  various  axial  rotations.  A  thrust-to-weight  ratio  of 
T/W  =  0.0125  is  initially  assumed;  the  effects  of  varying  T/W  are  also 
investigated. 

The  nature  of  the  physical  problem  itself  and  the  structure  of  the  resulting 
TPBVP  have  led  to  a  set  of  three  Cases  which  are  described  in  Table  1. 
Reference  will  be  made  extensively  throughout  the  rest  of  the  paper  to  Cases  1, 
II,  and  III.  Case  I  consists  of  transfers  between  fixed  points  on  the  initial  and 
terminal  orbits  whose  locations  are  specified  in  advance.  The  solution  of  the 
TPBVP  yields  the  optimum  thrust  orientation  history  as  well  as  the  transfer 
trajectory  and  fuel  expenditure.  Case  II  consists  of  transfers  from  a  speci¬ 
fied  point  on  the  initial  orbit  to  a  point  (not  specified  in  advance)  on  the  terminal 
orbit,  the  location  of  which  is  an  additional  part  of  the  problem.  Case  III 
consists  of  fuel  optimal  transfers  between  orbits  in  which  neither  the  departure 
nor  arrival  point  is  fixed  in  advance,  but  must  be  determined  as  part  of  the 
solution  to  the  problem. 
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NOTATION 

=  terminal  ellipse  eccentricity 
hj  =  terminal  ellipse  angular  momentum 
m  "  spacecraft  mass 
m  =  mass  rate  (constant) 
p^.  =  terminal  ellipse  semilatus  rectum 

r  =  radial  distance  to  spacecraft 
T  =  thrust  magnitude  (constant) 
tj  =  final  time 
t  =  time 


V  =  circumferential  component 

of  velocity 

V  =  radial  component  of  velocity 

V  =  total  velocity  V  =  +  Vy  j 

p,  =  gravitational  constant 

v  =  angular  position  of  spacecraft 
relative  to  perigee  of  terminal 
orbit 

<p  =  thrust  angle 

subscript  0  refers  to  initial  conditions 


Figure  1.  Geometry  and  Nomenclature  of  the  Orbit  Transfer 
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Table  1.  Orbit  Transfer  Cases 


Case 

Physical  Interpretation 

TPBYP  Interpretation 

1 

Fuel  optimal  transfer  from  a  fixed 
point  in  the  initial  orbit  to  a  fixed 
point  in  the  terminal  orbit 

Fixed  endpoint  optimization 
problem 

II 

Fuel  optimal  transfer  from  a  fixed 
point  in  the  initial  orbit  to  a  point 
(not  fixed)  in  the  terminal  orbit 

One  endpoint  fixed,  one  endpoint 
free  optimization  problem 

III 

Fuel  optimal  transfer  between  two 
orbits  in  which  neither  departure 
nor  arrival  point  is  fixed  in  ad¬ 
vance 

Two  free  endpoints  optimization 
problem 
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II.  FORMULATION  OF  THE  PROBLEM 


A.  DYNAMIC  EQUATIONS  AND  SPECIFIED  BOUNDARY  CONDITIONS 

Using  the  notation  defined  in  Fig.  1,  the  equations  of  motion  for  the  non- 
atmospheric  spherical  earth  gravitational  model  are: 


•  _  » 
r 


v2  „  . 

v  _  V  .¥•  ,  ■  T  sin  0 
r  r  Z  (m0-  |m|t) 

V  V 

v  -  r  v  ,  T  cos  0 
v  “  "  r  (mQ-  |m|t) 


(1) 


The  initial  conditions  are  specified  by: 

Pi 

r  (t  )  = - 1 _ 

'  0;  1  +  e.  cos  v(t_) 


*%)  =  VQspec  (Cases  1  and  11  onlY) 


Vr(t0)  = 


e.h.  sin  v(tQ) 


V‘0> 


(2) 


with  the  terminal  conditions  given  as 


r(tf) 

v(t£) 

Vr<V 


It  _ 

1  +  ef  cos  v(tf) 

■  Vfspec  (Case  1  only) 
efh£  sin  v(t£) 

=  Pf 


W 


(3) 


where  the  subscripts  i  and  f  refer  to  the  initial  and  final  orbits;  e^,  a^,  e^,  and 
a^  are  specified  orbital  parameters  (eccentricity  and  semimajor  axis);  and  h^, 
p..,  h^,  p^  are  the  corresponding  angular  momentum  and  semilatus  rectum  for 
each  ellipse. 

The  orbit  transfer  is  to  be  made  by  orienting  the  thrust  vector  at  the 
best  instantaneous  angle  <f>(t),  measured  clockwise  from  the  circumferential 
direction,  to  minimize  the  total  fuel  expenditure  for  the  maneuver.  Since  a 
constant  mass  flow  rate  is  assumed,  a  minimum  time  transfer  is  equivalent 
to  minimizing  the  fuel.  This  can  be  viewed  as  a  classical  Mayer  Problem  in 
the  Calculus  of  Variations,  and  the  necessary  conditions  for  optimality  are 
reviewed  in  the  next  section. 

B.  NECESSARY  CONDITIONS  FOR  OPTIMALITY: 

THE  MAYER  PROBLEM 

It  is  desired  to  minimize  the  criteria  function  of  the  form 

J  =  $(x,  t)|._.  (4) 

1  f 
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subject  to  the  differential  equation  constraints 


x  =  f_(x,  u)  ,  Xq  =  x(tQ)  15 ' 

where  x  and_f  are  n-vectors  and  u  is  a  scalar  function  of  time,  and  to  the  speci 
fied  terminal  constraints. 


M(x,t)  |t=t^  =  0  (6) 

where  M  is  an  m- vector  with  m  <  n.  If,  for  example,  the  complete  state 
vectoi  is  specified  at  the  terminal  point  (as  in  the  Case  I  orbital  transfer), 
then 


M 

M 

M 

M 


1 

2 

3 

4 


=  v,  - 


Pf 

i  +  e,  cos  v, 

f  fspec 

Vfspec 

efhf  Sin  Vfspec 

Pf 


(7) 


Consider  further  the  free  time  problem  where  the  final  time  tf  is  not  specified 
The  augmented  cost  function  is  constructed  as 


=  v(xf,  t£) 


+  n 


M(xf,tf) 


,  A  T 

A 


\T(t)[*  -i(x.u)]dt 


(8) 
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where  \(t)  is  an  n-vector  of  Lagrange  multiplier  functions  and  ii  is  a  constant 
m- vector  of  Lagrange  multipliers. 

Requiring  that  the  first  variation  of  J  vanish  yields  the  Euler- Lagrange 
equations 

-  -  \T  [costate  equations]  (9) 

and 

XT  =  ®  [optimality  condition]  (10) 


and  boundary  conditions 


\T(tf) 


8$  ,  „  T  9M 
8x  “  8x 


[‘•ransversality  conditions] 


(ID 


(XTx)  | 

l£ 


9$ 

at 


-  8t 


(12) 


In  principle,  Eq.  (10)  may  be  solved  at  each  time  point  for  the  instantaneous 
control 


u  =  u(x,\)  (13) 

and  substitution  may  be  made  into  Eq.  (9)  to  give  the  reduced  adjoint  equation 

h  =  z.(x>h)  (14) 


where 
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If  an  iterative  numerical  procedure  is  required  to  calculate  u(t),  starting 

T  2  2 

estimates  must  be  made  in  a  region  where  X.  3  _£/9u  <0,  as  required  by  the 
Weierstrass  Condition.  An  analytic  solution  must  satisfy  this  condition  as 
well. 

The  boundary  conditions  for  Eq.  (14)  are  given  at  the  final  time  by 
Eq.  (11),  whereas  the  initial  conditions  of  the  state  equations,  Eq.  (5),  are 
known.  Hence,  the  original  optimization  problem  of  Mayer  has  been  trans¬ 
formed  into  a  two-point  boundary  value  problem  (TPBVP).  If  Eq.  (5)  is  a 
set  of  non-linear  dynamic  equations,  an  analytic  solution  of  the  TPBVP  is  not 
possible  except  under  extremely  unusual  circumstances.  To  sum  up,  the 
unknown  quantities  are: 

no.  of  elements 

-0  =  -^V  n 

m 
1 

Total  n  +  m  +  1  unknowns 


n 

t. 


The  conditions  to  be  satisfied  are: 

no.  of  elements 


M 

terminal  constraints 

Eq.  (6) 

x(tf) 

transversality 

Eq.  (11) 

(XTx)| 

Eq.  (12) 

f 


Total  m  +  n  +  1  conditions 


Thus  it  is  seen  n  +  m  +  1  unknowns  are  available  to  satisfy  an  equal  number 
of  conditions;  the  system  is  determined. 

C.  REDUCTION  OF  THE  TRANS  VERSA  LIT  Y  CONDITIONS 

Since  the  minimum  fuel  orbit  transfer  using  a  continuous -burning 
constant  mass  flow  rate  engine  requires  a  minimum  time  traversal,  the 
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appropriate  Mayer  problem  cost  function  is  $(x(t),t)|,  =  t,. 

tf 

The  state  variables  may  be  assigned  as: 


x 


1 


=  r 


x 


2 


-  v 


X 


3 


=  V 


r 


x 


4 


The  specified  terminal  conditions  (Eq.  6)  for  Case  I  become 


-  x 

M2  =x 

M3  =  * 
M4  =  x 


Pf 


if  i  +  e^  cos  x 

2f  "  vfspec 

e^hj.  sin  x2£ 


3f 


4f 


2f 


(15) 


(16) 


Thus 


9M 

9x 


“f 

2 

Cif 


pfef  sin  x2f 


(l  +  e£  cos  x2£)‘ 
i 


cos  x 


2f 


Pf 


(17) 


and  since  9^/9x  =  (0,0,  0,0),  Eq.  (11)  gives 


Xl£=r|l  -  I 


p^e^  sm  x2£ 

1  +  e^  cos  x2£ 


12  '2 


X2f  =  ^2 


e£h£  cos  x 


2f 


k3f 


^2+tl3 


X4f  =  2  ^1  +  ^4 
Xif 


(18) 


The  Lagrange  multipliers  ri^ ,  r|2>  n3>  n4  have  no  specified  values  and  thus 
the  \j£,  X-2£,  X-3£,  X^£  are  undetermined  for  the  Case  I  boundary  conditions. 
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For  Case  II  transfers,  the  specified  boundary  conditions  are 


M .  =  x. 


1  If  1  f  cos 


M2  =  x3f 


efhf  sin  x2f 


=  x. , 

3  4f  x 


if 


(19) 


Thus 


Then  since  =  (0,  0,  0,  0), 


9Mt 

9Mt 

9M 

9M 

9xi 

Sx2 

9x3 

9x4 

9M2 

8M2 

(M 

2 

CO 

9M 

9x^ 

9x2 

9x3 

X4 

co 

1 

9M 

3 

9M3 

9M. 

9xi 

9x2 

9x3 

9x4 

x 


if 


p^ef  sm  x 


2f 


(i  +  ef  cos  x2iy 


e  hr  ros  x 

I  X 


2f 


(20) 


Xlf  =  ^1  +  2  ^3 
Xif 


Pfef  sin  x2f 


L2f 


(l  +  ef  cos  x2f) 


2  \ 


^f 


"N 


cos  X 


2f  _  V  Transver  sality  con- 
'  ditions  in  terms  of 
Lagrange  multipliers  t) 


X3£  ^2 


X4f  =  ^3 


J 


(21) 
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Thus,  the  terminal  coupling  of  the  state  and  adjoint  variables  can  be  expressed 
as 


\2£  + 


p^e^  sin  x^^ 

(i  +  ef  cos  x2f) 


COS  X, 


(22) 


If  the  transfer  is  constrained  to  depart  from  the  initial  ellipse,  but  at 
an  unspecified  exit  point,  an  additional  set  of  transversality  conditions  arises 
leading  to 


V20  + 


p.e. 

1 


sin  x 


20 


(l  +  e.  cos  *20) 


X 


e.h.  cos  x0n 

-Z- X4o)  *  '  ~  *  p. . .  So  '  ° 

10  '  1 


(23) 


Thus  Eq.  (22)  and  (23)  serve  to  couple  both  the  initial  and  final  state  and 
costate  variables,  and,  although  the  unknowns  and  cannot  be  determined 
directly  from  these,  they  provide  the  additional  information  necessary  to 
establish  the  optimal  values. 

D.  SCALING  OF  THE  VARIABLES 

The  state  variables  are  scaled  into  characteristic  units  in  the  following 
manner: 

(a)  Time:  t  =  (jx/R3)1  /2t  =  Ct  (24) 

where 

t  =  unsealed  time,  in  seconds 
[i  =  GM 

R  =  radius  of  the  earth 
G  =  universal  gravitational  constant 
M  =  mass  of  the  earth/satellite  system; 
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the  scaling  constant  C  has  radians /second  as  the  units  and  is  the  mean 
angular  velocity  of  a  surface  circular  satellite. 

(b)  Length:  r  =  r/R  (25) 

where  r  =  radial  distance  to  satellite  measured  from  the  center  of  the  earth, 
in  feet. 

(c)  Velocity:  V  =  V/RC  (26) 

1  /2 

where  V  =  velocity  of  satellite  in  feet/second  and  the  product  RC  =  (|J./R) 

and  is  thus  the  velocity  of  a  surface  circular  satellite  measured  in  feet/second. 

(d)  Time  normalization  (Long's  transformation): 

For  computational  purposes,  it  is  convenient  to  introduce  a  new  independent 
variable  s  by  the  transformation 


t  =  as 


(27) 


where  0  ^  s  £  1  and  a  is  a  constant. 

Thus  when  s  =  1,  a  =  t^  (the  scaled  final  time).  Integration  is  performed 
with  respect  to  s  for  a  fixed  number  of  steps,  since  its  range  is  constrained 
to  the  unit  interval.  For  free  (variable)  time  problems,  the  value  of  the 
unknown  parameter  a  is  determined  as  part  of  the  optimization  process.  The 
variable  s  is  referred  to  as  "normalized  time."  From  Eqs.  (24)  and  (27),  the 
differential  relationship  is  seen  to  be 


dt 


(28) 


and  the  state  and  costate  differential  equations  are  transformed  accordingly. 
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Corresponding  to  the  scaled  state  dynamic  equations  is  the  adjoint 
system 


V  X,  V2CX, 
v  2  +  v  3 


‘i  _2d 
r  R 


_2 

r 


2CX,  V  V  CX. 
3  r  v  4 

_3  '  -2 

r  r 


x2  =0 


X,  V  X, 
X, 

3  C  r 


X4  = 


2V  X_ 
v  3 


RCr 


V  K 

+  -4-i 

r 


(29) 


where  (  )'  =  d(  )/ds.  The  adjoint  variables  can  thus  be  nondimensionalized 
by  selecting  the  following  scaling  relations 


X 


i 


x2=x2/r 


(30) 


E.  REDUCTION  TO  STANDARD  FORM  FOR  COMPUTATION 

The  necessary  conditions  which  must  be  satisfied  by  the  optimal  orbital 
transfer  trajectory  constitute  a  nonlinear  two  point  boundary  value  problem 
which  in  general  does  not  admit  an  analytic  solution.  The  unknown  parameters 
which  must  be  determined  to  solve  the  TPBVP  are  the  initial  condition  vector 
for  the  costate  differential  system  X^  and  the  free  (unspecified)  final  transfer 
time  t^.  These  must  be  adjusted  to  satisfy  the  transversality  conditions  at  the 
initial  and  final  point,  when  applicable,  and  the  specified  elements  of  the 
terminal  state  vector. 
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For  computational  purposes,  it  is  convenient  to  define  a  vector  function 
r  in  the  following  manner: 


M 

------  (31) 

Ti 

T 

n-m+1 

where  M  represents  the  specified  conditions  on  the  state  vector  and  T.  consists 
of  the  transversality  conditions  (in  reduced  form).  Thus  since  the  functions 
in  Tare  evaluated  at  times  t^.  and  t^  (for  Case  III), 

iT  _  h.Q’%. £*  — f ’  ^ 

=  it^o’^-o’^f^-O’-O’  V’-f^-O’-O’  V*fcf^ 

=  Kxi(to)j^o'V  {32) 

where  x.(t  )  designates  the  unspecified  elements  of  =  2E(tg)*  A  vector 
of  unknowns  is  next  defined  as 


Z  = 


and  hence  the  problem  is  transformed  into  the  following: 

* 

find  Y.  such  that 


Xi(t0") 


(33) 


Kz*)  =  £ 


(34) 
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where  0  is  the  zero  vector  of  proper  dimension  and  is  the  optimal  value 
of  £  or  in  the  equivalent  scalar  form,  find 

min||r(£)|| 

I 

where  it  is  known  from  the  definition  of  Tthat  the  minimum  norm  is  zero. 

The  simplicity  of  the  relation  Eq.  (34)  is  deceptive  in  that  the  solution  y_ 
may  be  extremely  difficult  to  obtain.  There  are  several  reasons  for  this: 

(1)  since  y;  represents  initial  conditions  on  differential  equations,  the  value 
of  the  function  T  is  very  sensitive  to  small  changes  in  its  argument;  (2)  T  can 
be  evaluated  only  after  integration  of  the  state  and  costate  differential  equations 
from  tg  to  tj.,  which  involves  many  computations  and  numerical  integration 
errors  that  may  become  critical  here;  (3)  the  contours  of  |  |r|  |  are  almost 
always  highly  irregular  and  contain  narrow  channels  and  multiple  extrema. 
Nevertheless,  the  above  problem  may,  in  principle,  be  solved  by  numerical 
techniques  which  are  applicable  for  multivariable  function  minimization. 


17 


III.  DESCRIPTION  of  the  algorithms 


Two  finite  dimensional  optimization  techniques  which  were  employed  for 
solving  the  nonlinear  boundary  value  problems  for  the  optimal  orbit  transfer 
are  reviewed  in  this  section.  The  basic  principles  of  the  algorithms  will  be 
discussed  briefly.  For  further  details  the  reader  should  consult  the  references 

A.  MULTIPTF  SUBSTITUTION  POLYNOMIALS 

This  is  an  iterative  gradient- free  optimization  technique  which  acts  to 
generate  a  sequence  of  substitute  contour  systems  for  the  vector  function 
r(y).  Each  of  these  in  turn  is  treated  by  standard  techniques  applicable  for 

solving  simultaneous  nonlinear  algebraic  equations.  R.  F.  Jaggers  (Ref.  7) 
has  developed  the  following  method  for  calculating  and  applying  multivariable 

substitution  polynomials  in  dynamic  optimization  problems. 

The  basic  plan  is  to  represent  the  vector  function  r(y)  by  a  system  of 

second-order  multivariable  polynomials  as  follows: 


where  jr  =  (y. ,  ...  ,  is  the  vector  of  unknowns  and  the  a^.  are  constant 

coefficients  tc  'e  determined.  By  definition,  y^  =  1  so  both  constant  and 


Blank 
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first  order  terms  are  also  present  in  the  representation.  It  is  convenient  to 
define  a  transformation 


z.  =  (y.  -  E.)/A.  i  =  .  .  .  ,  n  (36) 

1  1  1  1 

where  E.  represents  the  best  estimate  of  y.‘  and  the  A.  are  selected  scaling 
i  1  ii 

parameters.  Thus  Eq.  (35)  becomes 


r 


1 


n  n 

=E2>! 

i=0j=i  J 


(i) 


z.z. 
i  j 


z.z. 
1  J 


n  iJ  1  J 

i=0  j=i 


(37) 


where  the  constants  a|^  are  related  to  the  ajk^,  but  only  the  former  need  to  be 
determined.  These  coefficients  may  be  calculated  without  recourse  to  a  matrix 
inversion,  using  the  algorithm  developed  by  Jaggers.  See  Ref.  7  or  8  for 
details  of  this  procedure. 

The  next  step  is  to  set  the  right  side  of  Eq.  (37)  to  the  zero  vector  and 

T 

find  the  vector  z  =  (zj,  ...  ,  z  )  which  satisfies  T(z)  =  0.  Denoting  by 
Z<k>  =  (z<k).  . .  ,  z^)T  the  solution  of  the  ktk  substitute  contour  system  (i.e. , 
the  km  iteration),  the  corresponding  value  for  y  is 


r!k*  =  A.z!k^  +  E. 


i  i 


i  =  1, 


n 


(38) 
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It  is  not  presumed  that  this  y^k)  is  the  exact  solution  to  the  original 
problem  since  it  can  be  only  as  precise  as  the  approximation  used  in  Eq.  (37). 
Instead,  is  treated  as  the  new  best  estimate  of  the  true  optimum  y  and, 

accordingly,  E.  =  y!k^  is  substituted.  New  polynomial  representations  are 
generated,  centered  about  this  new  estimate,  and  the  procedure  continues  as 
before.  Iterations  are  continued  until  r(y(k))~  0  to  the  required  tolerance. 
The  procedure  is  schematically  illustrated  in  Fig.  2  where  an  approximation 
for  a  single  function  of  two  variables  is  given.  The  solution  of  the  substitute 
contour  system  after  the  itk  iteration  (denoted  by  E  )  can  be  found,  for 
example,  by  using  the  Newton- Raphson  algorithm.  The  gradient  which  is 
required  can  readily  be  calculated  analytically  from  the  quadratic  form  as 
follows: 


Hl 

3zk 


k  =  i ,  .  .  .  ,  n 
l  =  1,  .  .  .  ,  n 


(39) 


where 


2  if  k  =  j 
1  if  k  f  j 


and 


if  k  >  j 


The  polynomials  are  established  by  evaluating  actual  function  T  at  test  points 
of  the  form 


yi  =  1 


E.  +  A. 

l  l 

E. 

l 

E.  -  A. 

l  l 


l  —  1 ,  .  .  .  ,  n 


(40) 
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and  the  representation  is  made  via  multivariable  interpolation  formulas  which 
agree  with  the  actual  function  at  these  test  points  only.  The  scale  parameters 
A^  should  be  selected  relatively  large  when  beginning  the  iterations  and  reduced 
under  program  control,  as  the  solution  is  approached,  to  refine  the  region  of 
interest.  The  desired  situation  occurs  when  the  test  points  "straddle"  the 
solution,  i.e.  ,  have  most  of  the  I\  change  sign  when  the  +A^  and  -  A.  test 
points  are  used.  Convergence  on  subsequent  iterations  will  then  usually  pro¬ 
ceed  quite  rapidly. 

The  principal  features  of  the  multiple  substitution  polynomial  method  are: 
(1)  the  gradient  of  the  function  T  is  not  required;  (2)  the  general  n  dimensional 
algorithm  allows  easy  computer  application,  and  programming  effort  and 
computer  storage  requirements  are  minimal;  and  (3)  convergence  is  rapid  in 
the  neighborhood  of  the  solution.  This  technique  was  found  to  perform  well 
for  parametric  studies  of  the  fixed  endpoint  orbit  transfer  problem. 

The  computational  difficulties  which  arise  in  applying  MSP  are  related 
to  the  accuracy  of  the  contour  approximation  and  the  iterative  solution  of  the 
substitute  contour  system.  When  the  approximation  is  relatively  weak,  a 
minimum  norm  solution  of  the  latter  can  sometimes  be  used  effectively  when 
the  zero  norm  does  not  exist.  An  efficient  iterative  method  for  solving  the 
generated  simultaneous  nonlinear  algebraic  equations  is  essential  for  success 
since  loss  of  convergence  here  causes  a  breakdown  of  the  whole  procedure. 

The  Newton- Raphson  algorithm  with  an  automatic  convergence  monitor  was 
found  to  give  good  results. 

B.  MARQUARDT'S  ALGORITHM 

This  is  a  gradient-oriented  method  which  seeks  to  combine  the  principal 
features  of  the  Taylor  series  method  (i.e.  ,  Newton- Raphson)  and  the  gradient 
method  (steepest  descent)  by  performing  an  adaptive  interpolation  between 
them.  The  basic  idea  is  to  solve  the  desired  system  r(y  )  =  £  by  generating 
a  correction  vector  Ay^  =  -  y^  which  minimizes  the  scalar  function 

S(Ay)  =  (r(y)  +  AAy)T(r(y)  +  A  Ay)  +  \(AyTAy  -  R2)  (41) 
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where  A  =  9r/9y  is  the  gradient  matrix  evaluated  at  the  point  y,  and  X  is  a 

scalar  Lagrange  multiplier  associated  with  a  correction  size  constraint 
T  2 

Ay  Ay  =  R  .  R  is  a  positive  constant  representing  the  radius  of  a  hyper¬ 
sphere  over  which  the  minimization  of  S  is  to  be  ^performed,  The  appropriate 
minimizing  correction  vector  corresponding  to  Eq.  (41)  is  readily  shown 
to  be 


Ay  =  -  (ATA  t  \I)_1ATr(y;)  (42) 

It  is  of  interest  to  observe  that  the  correction  step  given  by  the  Newton- 
Raphson  algorithm  is 


AZnr'-A'1^  (43) 

which  corresponds  to  the  limiting  case  X  =  0  in  Eq.  (42).  Furthermore,  the 

2,  T 

steepest  descent  direction  for  the  function  Fg  =  |  |r(y;)  |  |  =  r(y)  r(y)  is 

AySD  =  -  KATr(y)  (44) 

where  K  is  a  positive  constant.  This  is  approached  by  Eq.  (42)  as  X  -*•  °o. 

2 

Marquardt  (See  Ref.  9)  shows  that  (a)  |  1  Ay (M  |  |  is  a  continuously  decreasing 

2 

function  of  X  such  that,  as  X  -*•  oo,  |  j  Ay(X)  |  |  -*  0  and  (b)  the  angle  given  by 

-i  I  ^t^sd  \ 

v  =  cos  Irunr 

is  a  continuous  monotone  decreasing  function  of  X  such  that,  as  X  -*•  co,  y  -*•  0. 

It  follows  that,  since  Ay^^  is  independent  of  X,  the  vector  Ay  rotates  toward 
Ay^p  as  X  -*  w  .  In  constructing  the  algorithm,  the  strategy  for  adjusting  the 
parameter  X  is  as  follows: 
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At  the  k^k  iteration,  X^  must  be  chosen  such  that 
(k+i)  (k) 

r  far  )  <  I"  hr  ).  This  can  always  be  done  in  principle  if 

o  S 

the  steepest  descent  direction  is  defined  and  is  non-zero. 

(k) 

However,  this  may  necessitate  choosing  a  large  value  of  X.'  . 

(k) 

Xv  7  should  be  reduced  toward  zero  as  rapidly  as  possible  when 
conditions  exist  which  allow  convergence  with  the  Newton-Raphson 
correction  vectors.  Marquardt  defines  the  following  adaptive 
scheme  for  carrying  this  out: 

(a)  Let  v  >  1 


(b) 

(c) 


(k- 1 ) 

Let  X'  7  denote  the  value  of  X  from  the  previous  iteration. 
Initially  let  X^  =  i0-7%  say. 

Compute  Ts(£(k))  |  (k_t  j  and  rs(y:(k))  |  (k>1)  where 

r  (y(l))U  =  r  +  ^X^V))  and  Ay_  is  calculated  from 

S  o 


Eq.  (42).  The  superscripts  refer  to  the  iteration  number. 


(d)  Then  if  r  (y(k))  |  ,  ..  <  r  fc(k“1)),  let  X{k)  =X(k"i}/v. 

s  XlK“U/v  s 

If  rsfc<k)>lx(k-i)/v>rs(I(k-1))  and  rs(yW)|^(k.1)S  ryk-‘>), 
letx(k)  =X(k_1).  Ifrs(Z(k))|  {k_1}  >  rs(y(k_i))  and 


/v 


(k- 1 ) 

increase  \  by  successive  multiplications  of  v  until  for 


some  smallest  co, 


r.t«)l 


/i  4  \  -  T  (v(k  ^)  and  let 

^(k-1 )  j^ui  s  J‘ 


x(k)  _x(k-l)vo)^ 


Since  gradient  methods  are  sea’  i -dependent,  Marquardt  elects  to  scale 

T 

the  parameter  space  by  introduc’ag  the  transformation  matrix  D  =  diag  (A  A). 
The  scaled  correction  vector  corresponding  to  Eq,  (42)  then  becomes 


Ay  =  -  D“1/2(D_i  ^2ATAD“1/2  +  \I)-1D“  i/2ATr(y;) 


(45) 


and  this,  in  conjunction  with  the  procedure  above  for  manipulating  X.,  consti¬ 
tutes  the  complete  algorithm. 

The  principal  advantage  of  Marquardt1  s  method  lies  in  its  ability  to 
combine  the  features  of  steepest  descent  which  exhibit  good  starting  charac¬ 
teristics  from  an  initial  guess  and  the  Newton- Raphs on  procedure  which  gives 
fast  terminal  convergence  in  the  neighborhood  of  the  solution.  The  adaptive 

adjustment  cf  the  interpolating  parameter  X.  assures  that  the  procedure  will 

2 

not  diverge  and  a  finite  reduction  of  r  =  |  jr(y)  |  J  is  realized  on  each 
iteration. 

One  of  the  numerical-problems  which  has  been  fouiid  to  occur  is  a  very 

slow  rate  of  convergence  so  that  a  large  number  of  iterations  gives  only  a 

small  reduction  of  the  function  T  .  This  situation  occurs  when  the  calculated 

s 

X.  becomes  reduced  so  that  the  angle  between  Ay  and  Ay^  approaches  90  deg. 

Thus  very  little  improvement  per  iteration  is  obtained.  Furthermore,  a 

matrix  inversion  and  a  computation  of  the  function  gradient  is  required  for 

each  iteration.  This  latter  is  a  time-consuming  operation  for  dynamic 

optimization  problems  which  require  integration  of  differential  equations  for 

each  evaluation  of  the  fimction  F  .  Hence,  one  desires  maximum  efficient  use 

s 

from  each  gradient  computed.  Nevertheless,  Marquardt' s  method  is  a 
powerful  optimization  procedure  which  can  be  used  effectively  to  solve  two- 
point  boundary  value  problems.  Its  performance  compares  favorably  to 
most  f:rst-order  optimization  techniques  which  do  not  exploit  known  contour 
features  for  a  given  problem. 
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IV.  DISCUSSION  AND  RESULTS 


The  true  minimum  fuel  transfer  between  two  given  orbits  is  the  Case  III 
problem  in  which  transversality  conditions  are  derived  which  provide  addi¬ 
tional  information  for  determining  the  optimal  departure  point  from  the  initial 
orbit  and  the.  best  final  arrival  point.  Because  of  the  numerical  complexities 
involved,  it  is  convenient  to  treat  this  problem  as  a  sequence  of  subproblems 
with  the  boundary  conditions  completely  or  partially  specified,  i.e.  ,  the 
Case  land  Case  II  problems.  This  approach  has  several  advantages: 

(1)  computational  problems  are  reduced;  (2)  physical  insight  into  the  geometry 
of  low  thrust  transfers  is  provided  by  the  optimal  solution  to  the  subproblems; 
(3)  sensitivity  of  the  cost  function  to  non-optimal  departure  and  arrival  points 
is  available;  and  (4)  existence  of  multiple  extrema  can  be  demonstrated. 

A.  THE  LIMITING  CASE  I  PROBLEM:  PHYSICAL  BARRIER 

Consider  the  Case  I  transfer  starting  from  perigee  of  the  initial  orbit 
and  proceeding  to  a  specified  endpoint  on  a  given  co-apsidal  final  orbit  (see 
Fig.  3).  The  Multiple  Substitution  Polynomial  method  was  used  to  provide 
the  numerical  solution  for  the  minimum  fuel,  continuous  low  thrust  transfers 
to  various  points  on  the  terminal  orbit,  as  specified  by  the  final  true  anomaly 
Vj..  Figure  4  shows  the  time  variation  of  the  true  anomaly  on  the  optimal  tra¬ 
jectory  for  specified  of  270,  240,  210,  and  180  deg.  Time  is  normalized 
for  each  transfer  trajectory  by  dividing  by  the  total  transfer  time  required. 
Thus  each  trajectory  has  its  own  normalization  factor,  and  each  is  plotted 
against  the  same  fraction  of  its  total  time.  The  optimal  solutions  show  that 
the  true  anomaly  is  nearly  a  linear  function  of  time  for  these  low- eccentricity 
orbits,  and  thus  a  minimum  time  transfer  acts  to  minimize  the  total  change 
in  true  anomaly. 

The  Case  I  problem  specifies  the  Av  =  v ^  for  the  transfer  and  it  is 
evident  that  there  exists  some  which  represents  a  physical  limit  of  the 

thrust  capability  of  the  vehicle.  For  a  fixed  v^,  no  Case  I  solution  exists  for 

Vf  <  without  a  large  jump  in  the  cost  function  caused  by  a  transfer  of 

min 
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more  than  one  revolution.  Thus,  the  Case  II  solution  represents  the  limit 
of  the  Case  I  fixed-point  transfers  as  approaches  v'fm£n*rom  above.  The 
Vfmin  *S  a  Physical  constraint  such  that  no  sample  points  can  be  evaluated 
during  the  iteration  process  which  requires  <  v£m^n* 

The  sequence  of  substitute  contours  for  r(y)  generated  by  the  MSP 
method  was  found  to  provide  rapid  convergence  for  the  specified  final  points 
taken  over  the  range  180  deg  -  -  330  deg,  with  fixed  at  perigee.  Approxi¬ 

mately  5  iterations  were  required  for  complete  convergence  for  a  30  deg 
change  of  final  true  anomaly.  However,  as  was  reduced  further,  the  region 
of  convergence  diminished  sharply  and  required  a  greater  number  of  iterations 
to  produce  neighboring  solutions  for  1  to  5  deg  changes  in  v^.  The  method 
became  completely  ineffective  for  £  158  deg.  The  actual  contours  of  r(y) 
evidently  could  not  be  approximated  sufficiently  accurately  because  of  the 
effect  of  the  physical  border. 

Marquardt's  method  was  proved  to  be  useful  for  extending  the  solution 
region  for  v^<  158  deg.  This  technique  acts  to  minimize  the  sum  of  squares 
of  the  components  of  T  and,  through  its  adaptive  manipulation  of  the  step  size 
parameter,  local  non-zero  minima  were  sometimes  encountered.  These  had 
no  apparent  physical  significance  and  often  arrested  the  convergence  process. 

An  increase  in  numerical  sensitivity  was  observed  as  the  border  was  approached, 
and  the  trajectory  integration  range  partitioning  was  increased  from  100  to 
200  steps  to  preserve  numerical  accuracy.  The  rates  of  change  of  the  optimal 
initial  conditions  X^q,  ^3q»  ^40  were  seen  to  experience  a  rapid  rate  of  increase 
as  was  decreased  toward  the  minimum  value  (See  Fig.  5).  This  caused  a 
reduction  in  the  region  of  convergence  for  neighboring  solutions  as  v ^  was 
varied. 

Calculations  of  the  gradient  by  finite  differencing  was  found  to  require 
a  careful  adjustment  of  the  perturbation  parameter  6  especially  for  solutions 
near  =  155  deg.  Figure  6  illustrates  the  rapid  fluctuation  in  the 
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approximation  to  the  partial  derivative,  as  computed  from 


Al' 


+  6'X30'X40’a^  “  r^10'X30,X40’a^ 


AX. 


10 


Ar  _  r^10,?"30,X40,a  +  6^  "  r^10,X30,X40,a^ 
Aa  "  6 


(46) 
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with  6  ranging  from  10  to  10  .  The  curves  show  the  finite  difference 

approximations  evaluated  at  the  solution  for  =  155  deg  where,  in  this  case, 
the  scalar  r  is  defined  by 


4 

r=I>il  (47) 

i=i 

The  significant  features  to  note  in  Fig.  6  are  the  restricted  range  of  6  over 

which  the  partial  derivative  approximations  remain  relatively  constant  and 

-  8 

the  rapid  falloff  of  the  curves  for  6  5  10  .  This  latter  point  represents  the 

loss  of  numerical  accuracy  from  integration  truncation  and  roundoff  effects, 
and  it  fixes  the  lower  limit  for  5.  Comparison  with  Fig.  7  shows  the  relative 
behavior  of  the  derivative  approximations  taken  about  the  solution  points  for 

v,  =  240  and  100  deg.  Notice  that  the  range  over  which  the  curves  remain 

1  -7  -4 

constant  extends  through  10  £  6  £  10  ,a  considerable  increase  over  the 

7  6 

restricted  range  of  good  approximation  10  S  6  <  10  required  for  the  = 
155  deg  example.  Although  no  generalizations  concerning  this  phenomenon 
are  evident,  it  suggests  the  need  for  a  sensitivity  analysis  of  this  kind  when 
using  finite  difference  gradient  approximations  for  dynamic  optimization 
problems. 
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Since  the  absolute  value  definition  is  used  for  r  [see  Eq.  (47)],  the 

value  of  the  partial  derivatives,  as  seen  in  Figs.  6  and  7,  evaluated  at  the 

solution  are  not  zero  but  are  related  to  the  slopes  of  the  function  Tin  the 

direction  of  the  coordinate  axes.  As  is  taken  near  the  border  of  the  Case  I 

feasible  region,  these  slopes  approach  zero.  The  determination  of  the  true 

limiting  point  v ^  .  was  made  using  the  Case  II  problem  formulation  with  the 

single  transversality  condition.  The  optimal  Case  II  transfer  from  perigee 

was  found  to  arrive  at  the  final  orbit  with  v,  -  v,  .  =  90.  6  deg. 

f  f  mm  b 

B.  PARAMETRIC  CASE  II  SOLUTIONS 

Marquardt's  method  was  used  to  generate  a  sequence  of  Case  II  (final 
endpoint  free)  solutions,  with  the  specified  starting  points  taken  90  deg  apart. 
Typically,  7  iterations  were  required  for  convergence  to  the  new  solutions  when 
the  initial  takeoff  point  v.  was  changed  by  a  90-deg  interval.  Figure  8  shows 
the  optimal  transfer  trajectories  given  as  a  polar  plot  of  altitude  versus  true 
anomaly.  The  arrows  indicate  the  thrust  direction  and  the  cross  lines  mark 
the  point  at  which  the  thrust  changes  sign.  If  the  optimal  transfer  times 
(and  equivalently  fuel  expenditure)  are  plotted  as  a  function  of  the  specified 
initial  true  anomaly,  the  curve  of  Fig.  9  results.  This  shows  that  there  is  a 
unique  minimum-fuel  transfer  possible  between  these  given  co-apsidal  ellipses 
and  furthermore  indicates  the  magnitude  of  the  fuel  penalty  of  initiating  the 
orbit  transfer  at  other  than  the  optimal  point.  The  true  minimum  was 
determined  by  solving  the  Case  III  problem  with  the  two  transversality  condi¬ 
tions,  and  this  gave  =319  deg  as  the  best  departure  angle.  The  optimal 
transfer  trajectory  lies  in  the  mutual  perigee  region  and  is  indicated  in 
Fig.  8.  The  sequence  of  Case  II  solutions  served  to  isolate  the  proper  region 
for  the  Case  III  transfer  and  provided  sufficiently  close  initial  adjoint  condi¬ 
tions  so  that  the  five- dimensional  optimization  program  gave  rapid  convergence. 
The  Case  I  and  II  problems  required  only  a  four- dimensional  search  as 
formulated. 
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TRANSFER  TIME,  sec 


c. 


AXIAL  ROTATIONS  AND  THE  MULTIPLE  CASE  III  SOLUTIONS 


An  example  of  optimal  low-thrust  transfer  between  intersecti,  g  orbits 
was  studied,  using  the  same  orbits  employed  previously  but  rotated  with 
respect  to  each  other  by  Aiv  =  90  deg.  The  intersection  points  were  found  to 
occur  at  v  =  84  deg  and  155  deg  as  measured  from  perigee  of  the  initial 
ellipse. 

Parametric  solutions  of  the  Case  II  problem  were  conducted,  taking 
selected  starting  points  around  the  initial  orbit.  Severe  computational 
problems  arose  in  the  neighborhood  of  the  intersection  points  where  very  fast 
rates  of  change  of  the  optimal  adjoint  initial  conditions  were  observed.  The 
nature  of  the  computational  barriers  is  seen  in  Fig.  10.  The  solutions  in  the 
vicinity  of  \>.  -  130  deg  were  especially  difficult  to  obtain  since  the  region  of 
convergence  from  neighboring  solutions  was  quite  limited.  This  problem  was 
surmounted  by  increasing  the  numerical  integration  accuracy;  introducing 
more  precise  gradient  calculations  using  symmetric  differencing;  and 
reversing  the  direction  of  the  vc  parameter  variation.  In  the  regions  outside 
of  the  double  spike  area  shown  in  Fig.  10,  an  average  of  five  iterations  was 
sufficient  for  convergence  for  a  30  deg  change  in  the  specified  ic,  using 
Marquardt's  method. 

The  results  of  the  Case  II  parametric  study  for  the  90-deg-rotated 
geometry  are  displayed  in  Fig.  11  which  is  a  plot  of  the  non-dimensional 
transfer  time  t^  (scaled  by  a  constant  factor),  shown  as  a  function  of  specified 
departure  true  anomalies.  The  existence  of  multiple  extrema  is  clearly 
displayed  by  this  technique  of  using  parametric  studies  of  the  Case  II  problem. 
The  corresponding  optimal  adjoint  variable  initial  conditions  from  Fig.  10 
can  thus  be  selected  in  the  proper  regions  for  finding  the  multiple  minima 
for  the  Case  III  problem,  in  which  neither  endpoint  is  specified.  This  latter 
problem  requires  a  five-dimensional  search  procedure,  but  good  starting 
values  are  available  through  the  above  technique. 

Figure  12  shows  three  stationary  solutions  for  the  Case  III  orbit  transfer 
with  Aco  =  90  deg  and  displayed  as  a  polar  altitude  plot.  The  lesser  minimum 
occurs  for  the  transfer  trajectory  beginning  at  v  -  62  deg  on  the  initial  ellipse; 


Figure  10.  Case  II:  Initial  Values  of  Adjoint  Varia¬ 
bles  versus  True  Anomaly  of  Departure 
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Figure  12.  Case  III:  Polar  Plot  of  Altitude  versus  Departure 
Angle  for  90-deg  Rotated  Ellipses 


the  solution  starting  at  =  130  deg  is  a  local  minimum  only.  The  remaining 
trajectory  represents  a  relative  maximum.  A  second  relative  maximal 
extremum  is  located  between  the  minimizing  solutions  but  is  not  shown. 

D.  PROPERTIES  OF  THE  OPTIMAL  LOW  THRUST  TRAJECTORIES 

The  optimal  Case  III  continuous  low  thrust  orientation  histories  for  a 
variety  of  relative  orbit  axial  rotations  Aw  are  shown  in  Fig.  13.  The  thrust 
orientation  angle  <p  is  measured  clockwise  from  the  circumferential  direction, 
and  each  steering  history  is  plotted  with  respect  to  the  normalized  time  of 
its  own  transfer.  A  relatively  rapid  change  in  thrust  direction  is  seen  to 
occur  at  the  halfway  point  of  the  transfer  (T  =  0.  5)  for  the  longer  trajectories, 
as  shown  by  the  Aw  =  0  and  30  deg  curves.  From  empirical  results,  the 
optimal  transfer  trajectory  appears  symmetrical  about  the  point  of  closest 
approach  of  the  two  orbits.  The  rate  of  change  of  the  thrust  angle  decreases 
as  the  transfer  distance  decreases,  as  indicated  by  the  Aw  =  60  deg  example, 
although  symmetry  is  preserved.  The  orbits  intersect  for  Aw  -  84  deg  and 
multiple  extrema  were  found  to  exist  for  the  A^.  -  90  deg  case  examined.  The 
optimum  thrust  orientation  history  corresponding  to  the  least  fuel  transfer 
is  shown  for  the  90  deg  rotation  case. 

The  effect  of  the  thrust-to-weight  ratio  (T/W)  on  the  relative  fuel  cost 
and  transfer  time  is  shown  in  Fig.  14.  It  is  seen  that,  as  the  thrust  capa¬ 
bility  for  a  given  vehicle  is  reduced,  the  transfer  maneuver  becomes  more 
economical  from  the  fuel  cost  standpoint  but  at  the  expense  of  an  increase  in 
the  transfer  time.  The  curves  shown  were  computed  for  the  Case  III  problem 


THRUST  ANGLE  <f>, 


Figure  13.  Case  III:  Optimal  Thrust  Orientation 
versus  Time  for  Various  Aw 
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FUEL  EXPENDITURE,  PERCENT  OF 
INITIAL  VEHICLE  WEIGHT 


Figure  14.  Case  III:  Transfer  Time  and  Fuel  Expenditure 
versus  Thrust-to  Weight  Ratio 
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V.  SUMMARY  AND  CONCLUSIONS 


Some  computational  aspects  of  the  optimal  continuous  low  thrust  orbit 
transfer  problem  have  been  presented.  Numerical  solutions  were  obtained 
v;.a  the  method  of  Multiple  Substitution  Polynomials  and  Marquardt's  method, 
indicating  their  applicability  for  solving  optimal  control  problems  of  this 
type.  Formulation  of  the  problem  for  numerical  solution  was  given,  showing 
how  a  sequence  of  subproblems  can  be  used  to  indicate  the  existence  and 
location  of  multiple  extrema.  Computational  problems  relating  to  the  break¬ 
down  of  the  numerical  technique  on  restriction  of  the  region  of  convergence 
were  discussed. 

A  physical  barrier  which  restricts  sample  points  was  found  for  fixed 
endpoint  transfer  problems,  and  the  very  rapid  variation  in  the  optimal 
adjoint  initial  conditions  near  this  region  was  shown,  which  accounts  for  the 
severe  convergence  difficulties  encountered  here. 

Parametric  studies  of  the  Case  II  transfers  between  intersecting  orbits 
were  hampered  by  nearly  singular  behavior  in  the  Xq  vs  characteristic 
which  created  a  computational  barrier  for  isolating  the  multiple  solutions  for 
this  problem.  This  phenomenon  was  not  found  for  non-intersecting  orbits. 

Sensitivity  problems  for  calculating  the  function  gradient  by  finite 
differencing  techniques  were  encountered,  which  caused  loss  of  convergence 
of  Marquardt's  method  in  certain  transfer  regions.  The  critical  adjustment 
of  the  differencing  parameter  6  was  demonstrated  at  several  Case  II  solution 
points. 

The  optimal  continuous  low  thrust  orbit  transfer  trajectory  for  the 
free  endpoints  case  was  empirically  found  to  be  symmetrically  located  with 
respect  to  the  point  of  closest  approach  of  the  two  orbits.  A  rapid  thrust 
direction  reversal  was  seen  to  occur  at  a  time  equal  to  approximately  one 
half  of  the  ultimate  transfer  time  when  the  orbits  were  nearly  co~apsidal. 
This  corresponds  to  the  longer  duration  optimal  transfer  maneuvers  among 
the  cases  studied  and  the  form  of  this  solution  is  in  agreement,  for  example, 
with  the  continuous  low  thrust  Earth- to-Mars  rendezvous  results  reported  by 
Melbourne  and  Sauer  (Ref.  4)  and  Zimmerman  et  al  (Ref.  5). 
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From,  a  sequence  of  Case  III  solutions  with  various  T/W  ratios,  it  is 
shown  that  a  vehicle  with  a  larger  T/W  capability  requires  more  fuel  expendi¬ 
ture  in  transferring  between  two  given  orbits  than  does  one  with  a  smaller 
T/W.  The  total  time  for  the  transfer  maneuver  constitutes  the  major  tradeoff. 
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1  13  ABSTRACT 

The  purpose  of  this  paper  is  to  examine  the  computational  aspects  of 
the  minimum  fuel  continuous  low  thrust  orbit  transfer  problem  and  to 
display  characteristic  numerical  features  introduced  by  various  physical 
constraints.  Minimum -fuel  orbit  transfer  by  low  thrust  is  typical  of 
many  problems  in  optimal  control  which  result  in  a  two  point  boundary 
value  problem  which  must  be  solved  by  some  iterative  numerical  pro¬ 
cedure.  Two  techniques,  Multiple  Substitution  Poj.ynomin.als  (MSP)  and 
Marquardt's  method,  are  shown  to  be  applicable  to  this  task,  and  a 
detailed  analysis  is  made  of  the  behavior  of  these  methods  in  the  context 
of  the  low  thrust  problem.  A  variety  of  subproblems  is  considered  with 
parametric  variation  of  endpoints,  thrust-to -weight  ratio,  and  orbit 
axial  orientation.  A  physical  barrier  is  found  which  restricts  sample 
points  in  certain  limiting  case  fixed  endpoint  transfers.  The  existence 
of  multiple  stationary  solutions  is  shown  for  the  case  of  intersecting 
orbits,  and  the  nearly  singular  behavior  in  that  region  is  investigated. 
Numerical  results  for  several  transfers  are  found  to  compare  with  similar 
results  reported  elsewhere. 
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